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ABSTRACT 


Numerical  solutions  using  the  Hartree 
characteristic  method  for  the  propagation  of  spherical 
blast  waves  in  a  detonation  gas  are  presented  in  the 
present  paper.  Two  models  have  been  studied,  the  reacting 
blast  wave  (R-B-W)  model  and  the  finite  kinetic  rate  (F-K-R) 
model.  The  R-B-W  model  assumes  the  reaction  rates  to  be 
infinitely  fast,  and  the  reactive  shock  wave  is  treated  as 
a  discontinuity.  In  the  F-K-R  model  the  double  front 
approximation  is  used.  Empirical  relationships  for  the 
dependence  of  the  induction  time  on  reactant  concentrations 
and  thermodynamic  states  used  in  the  present  calculation 
are  taken  from  the  experimental  work  of  White  and  Strehlow. 
the  location  of  the  reaction  front  is  determined  using  the 
quaSi-steady  kinetic  approximation  as  proposed  by  Strehlow. 
For  the  R-B-W  model,  the  present  numerical  solution  is 
used  as  a  reference  to  assess  the  accuracies  and  the  domains 
of  validity  of  the  various  existing  analytical  solutions 
(i.e.,  quasi-similar  method  of  Oshima,  perturbation  method 
of  Melnikova  and  Sakurai  and  the  density  profile  method  of 
Pdrzel).  It  is  found  that  the  density  profile  method  of 
Pprzel  is  the  most  accurate  as  well  as  having  the  widest 
range  of  validity  for  both  the  parameters  at  the  shock  and 
the  flow  profiles  behind  the  blast.  The  present  solution 
also  confirms  the  result  of  Le.vin  and  Chernyi  in  that  the 
C-d  state  is  reached  at  a  finite  radius  for  cylindrical  and 


spherical  waves.  The  asymptotic  formula  of  Levin  and 
Chernyi  also  predicts  the  correct  value  for  the  C-J  radius. 
However,  it  is  found  from  the  present  solution  that  the 
J-T-Z  similarity  solution  as  assumed  by  Levin  and  Chernyi 
is  not  recovered  in  the  limit  as  the  C-J  state  is  reached. 
Instead  the  solution  tends  towards  the  planar  solution.  A 
new  asymptotic  decay  law  is  proposed  and  found  to  be 
independent  of  the  flow  gradients  to  the  first  order  in 
and  is  identical  to  Levin  and  Chernyi ‘s  formula.  Hence 
it  clarifies  the  correctness  of  Levin  and  Chernyi 's 
asymptotic  law  which  is  based  on  an  incorrect  assumption 
of  the  asymptotic  form  of  the  solution  since  the  C-J  radius 
is  independent  of  the  flow  gradient.  The  finite  kinetic 
Solution  reveals  all  the  essential  qualitative  features  of 
a  sphsrical  combustion  wave.  Quantitative  predictions  of 
thfc  sht-ck  and  reaction  front  trajectories  as  well-  as  the 
critical  ^.-iiergy  for  direct  initiation  are  correlated  with 
previous  experimental  data.  Based  on  a  blast  wave  model  for 
a  multiheaded  detonation  front,  the  present  F-K-R  solution 
is  used  to  predict  the  transverse  wave  spacings. 
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1 .  Introduction 

The  propagation  of  blast  waves,  in  a  detonating 
gas  has  received  considerable  attention  In  the  past 
decade  The  majority  of  the  theoretical  studies 

are  of  the  analytical  type  where  the  classical  point  blast 
similarity  solution  of  Taylor  and  Sedov  is  extended 

.  to  account  for  the  energy  release  due  to  chemical  reactions* 

In  these  reacting  blast  wave  (R-B-W)  models,  the  detonation 
front  Is  assumed  to  be  a  discontinuity  (I.e.,  Infinite  reaction 
rates)*  hence  all  the  non-similar  techniques  developed 

fOr  non-reacting  blast  waves  with  counter-pressure  effects  are 
directly  applicable  to  the  R-B-M  model  All 

these  non-similar  solutions  of  the  R-B-W  model  account  for 
small  departures  from  the  strong  point  blast  similarity 
solution  and  hence  break  down  as  the  initially  overdriven 
blast  wave  approaches  the  Chapman-Jouguet  (C-J)  condition. 

Their  validity  Is  restricted  to  the  regime  where  the  shock 
radius  Rg  Is  small  as  compared  to  the  characteristic  explosion 
length  R^  (i.e.,  Rj/Rq  <  !)•  It  is  one  of  the  aims  of  the 
present  paper  to  assess  the  extent  of  validity  of  the 
non-similar  analytical  solutions  In  terms  of  the  numerical 
solution  of  the  R-B-H  model  presented  herein.  The  asymptotic 
motion  as  the  blast  approaches  the  C-J  state  has  been  analysed 
by  Levin  and  Chernyi  with  the  rather  surprising  conclusion 
that  while  a  planar  reacting  blast  decays  to  a  C-J  wave 
asymptotically  (l.e^,  M  as  R_  -*■  »)  a  cylindrical  or 

spherical  reacting  blast  wave  approaches  the  C-J  state  at  a 
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finite  radius  (I.e.,  as  *»•  Rqj)*  Levin  and 

Chernyl's  asymptotic  decay  law  Is  based  on  the  assumption 
that  as  the  blast  wave  approaches  the  C-J  state  the  solution 
tends  to  the  classical  Jouguet  ^^®^-Taylor  ^^®^-Zeldov1ch 
(J-T-Z)  similarity  j^olutlon  for  a  diverging  C-J  detonation. 

It  Is  not  at  all  obvious  that  the  flow  gradients  behind  a 
decaying  blast  should  steepen  and  become  Infinite  at  the  C-J 
statR  as  described  by  the  J-T-Z  similarity  solution.  Hence 
In  th-si  present  paper  it  seemed  worthwhile  to  verify  this 
asymptotic  dec«y  law  of  Levin  and  Chernyi  by  comparing  it 
with  an  exact  rrar.^rlcal  solution  of  the  R-B-W  model.  It  Is 
particularly  Interesting  to  see  If  the  flow  gradients  behind 
the  decaying  blast  do  Indeed  steepen  and  take  on  the  J-T-Z 
profiles  as  the  C-J  state  Is  reached. 

In  an  effort  to  understand  the  chemical  reactions  In 
the  transient  hydrodynamic  flow  structure  of  a  decaying  blast 
wave,  theoretical  studies  on  a  finite  kinetic  rate  (F-K-R) 
model  have  been  made  These  studies  are  of  an 

analytical  nature  where  the  hydrodynamics  and  the  chemical 
kinetics  are  decoupled  by  assuming  that  the  heat  releasew  by 
chemical  reactions  has  a  negligible  Influence  on  the  motion 
of  the  blast  v^ave.  Based  on  this  assumption,  the  hydrodynamic 
flow  Structure  can  be  obtained  from  the  non-similarity  solution 
Of  a  non-reacting  blast,  and  with  the  flow  structure  known  the 
kinetic  rate  equations  can  then  be  Integrated  to  establish  the 
reaction  profiles  and  Induction  time  behind  the  blast  wave. 

Such  analytical  solutions  exclude  a  priori  the  all-important 
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coupling  mechanisms  between  the  energy  released  by  chemical 
reactions  and  the  hydrodynamic  flow  field.  Hence  these 
solutions  cannot  throw  light  on  the  interesting  experimental 
observations  such  as  the  re-establ  ishrr.ent  phenomenon  and 
the  existence  of  a  critical  energy  for  direct  initiation. 

The  present  paper  also  describes  the  result  of  a  numerical 
solution  of  the  F-K-R  model.  The  theoretical  predictions  of 
the  various  propagation  regimes  of  a  reacting  blast,  the 
dependency  of  the  critical  energy  for  direct  initiation  on 
the  initial  condition  of  the  explosive  system  as  well  as  the 
dimensions  of  transverse  wave  spacings  of  multiheaded  detonation 
waves  using  the  blast  wave  model  are  compared  with  existing 
experimental  results. 


2.  Theoretical  Models  and  Method  of  Solution 

The  basic  model  Is  that  of  a  point  spherical  blast 
wave  propagating  in  a  detonating  gas.  In  the  reacting  blast 
wave  (R-B-W)  formulation,  the  reaction  rates  are  assumed  to  be 
infinitely  fast  so  that  In  effect  the  chemical  reactions  release 
the  energy,  Q  per  unit  mass,  right  at  the  shock  front  which  is' 
treated  as  a  discontinuity.  In  the  fin  te  kinetic  rate  (F-K-R) 
formulation,  the  two-front  approximation  used  by  Bishimov, 
Korobeinikov  and  Levin  is  adopted.  The  two-front  model 

assumes  that  the  chemical  energy  Q  is  released  at  the  reaction 
front  only.  The  reaction  front  trails  behind  the  leading  shock 
front  by  an  induction  distance.  The  location  of  the  reaction 
front  is  determined  by  the  standard  quasi-steady  kinetic  method 


. . . . . . . . . . . . . . . . . . . . . .  . alilaai . . . . . . . . . wo . .  . . . I . . 


in  his  study  of  the 


first  used  by  Strehlow 

propagation  of  a  reaction  shock  wave  in  a  converging 

(22) 

channel  and  later  used  by  Bach,  Knystautas  and  Lee  '  , 

Lundstrom  and  Oppenheim  and  Bishimov,  Korobeinikov 
and  Levin  in  their  studies  of  the-  F-K-R  model.  The 

empirical  experimental  relationship  for  the  dependence  of 
the  induction  time  on  the  fuel -oxygen  concentrations  and 
the  thermodynamic  variables*  i.e., 

m  /-m  0 

(og^.  [O2]  [fae/]  )  =  A  +  (2.1) 

are  used.  In  the  numerical  solution  of  the  gasdynamic 
equations,  the  Hartree  characteristic  method  as  used  by 
Chou  and  Huang  in  their  numerical  solution  for  a 
non-reacting  blast  wave  is  adopted  in  the  present  work. 

The  perfect  gas  constant  y  assumption  is  used  throughout 
the  present  study.  The  details  of  the  basic  conservation 
equations,  the  Rankine-Hugoniot  relationships  for  the 
reacting  shock  front  as  well  as  the  method  of  computation 
using  the  Hartree  characteristic  scheme  have  already 
appeared  in  journal  publications.  It  suffices  here  to 
mention  the  appropriate  references 

3.  Comparison  with  Analytical  Solutions  for  the  R-B-W  Model 

Existing  analytical  solutions  for  the  R-B-H  model 
are  based  on  the  following  non-similar  techniques  developed 
to  account  for  counter-pressure  effects  in  non-reacting  blast 


wavesj  i)  the  quasi-similar  technique  of  Oshima,  ii)  the 
perturbation  solution  of  Melnikova  and  Sakurai  and  iii)  the 
power*"1aw  density  profile  method  of  Porzel .  Oshima 's  method 
was  first  used  by  Lee  for  the  R-B-M  model  in  1965  and 
later  by  Struck  and  Brossard  The  method  is  based 

on  the  assumption  that  the  flow  field  is  locally  similar. 

Hence  the  partial  differential  gasdynamic  equations  are 
reduced  to  ordinary  differential  equations  based  on  the 
instantanoous  shock  strength.  Oshima's  method  conserves 
the  total  energy,  hence  shock  trajectories  and  the  variation 
of  shock  strength  with  radius  are  adequately  described. 

However,  the  flow  structure  behind  the  blast  wave  is  poorly 
predicted  since  the  method  fails  to  conserve  the  total  mass 
enclosed  by  the  blast  at  any  instant.  The  method  of  Melnikova 
and  Sakurai  is  based  on  perturbing  the  similarity  solution 
of  Taylor  and  Sedov  for  the  strong  point  blast.  The  solution 
is  written  in  the  form  of  a  power  series  in  l/M^^.  This 
perturbation  method  was  used  by  Korobeinikov  and  Bach 
and  Lee  for  the  R-B-W  model.  The  method  is  rigorous  but 
due  to  the  asymptotic  nature  of  the  perturbation  expansion, 
the  Solution  diverges  rapidly  as  the  blast  decays  to  the  C-J 
condition.  The  power-law  density  profile  method  was  originally 
suggested  by  Porzel  and  developed  by  Rae  and  Bach 

and  Lee  The  method  assumes  the  form  Of  the  density 

profile  behind  the  blast  wave  to  be  given  by  a  simple  power^law 
form*  Then  from  the  conservation  of  mass  and  momentum  equations, 
the  forms  for  the  particle  velocity  and  the  pressure  profiles 
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can  be  obtained.  The  complete  solution  is  obtained  using 
the  known  profiles  in  the  energy  integral.  This  method  was 
applied  to  the  R-B-W  model  by  Bach  and  Lee  '  '  and  found  to 
be  an  extremely  accurate  method  even  as  the  blast  approaches 
the  condition.  In  comparing  the  present  solution  with 
these  various  analytical  solutions,  the  particular  case  of 
stoichiometric  C2*^2"^2  initial  pressure  of 

100  torr  and  at  room  temperature  is  used.  Thus  Fig.  1 
compares  the  variation  of  shock  strength  n'  = 
with  Shock  radius  Rg/Rg  fTom  the  various  analytical  solutions 
with  the  present  numerical  solution.  It  can  be  seen  that 
for  the  relatively  strong  reacting  blast  regime  (i.e., 

0  1  n'  <  .25)  all  the  analytical  solutions  are  almost 
identical  to  the  numerical  solution.  Departures  from  the 
exact  numerical  solution  occur  in  the  range  ,25  1  n'  11. 

The  density  profile  method  of  Porzel  yields  a  solution  that 
follows  closely  to  the  numerical  solution  throughout  the 
entire  range  0  1  n'  11. 

In  Fig.  2,  the  pressure  profiles  behind  the  reacting 
blast  wave  for  various  blast  radii  from  the  analytical 
solutions  are  compared  with  the  present  numerical  solution. 
For  the  relatively  strong  blast  regime  (i.e.^  0  1  n*  1  .4) 
the  profiles  from  the  perturbation  and  the  Porzel  solutions 
are  almost  identical  to  the  exact  numerical  solution.  For 
the  region  close  to  the  shock  front,  the  profiles  using 
0shima*s  method  are  identical  to  those  fron  the  numerical 
solution;  however,  they  begin  to  deviate  progressively  from 


it  as  one  approaches  the  center  of  symmetry.  For  the 
regime  of  moderate  shock  strength,  the  perturbation  solution 
retains  the  close  agreement  with  the  exact  solution.  Porzel's 
method  yields  a  good  description  of  the  profile  near  the 
front  but  predicts  a  lower  value  for  the  pressure  near  the 
center  of  symmetry.  All  the  analytical  solutions  break  down 
in  the  asymptotic  regime  as  except  the  method  of 

Porzel  which  still  yields  an  adequate  approximate  description. 

The  density  distribution  behind  the  decaying  reactive 

blast  is  shown  in  Fig.  3.  Except  for  the  very  strong  blast 

regime  where  n'  £  .25,  the  perturbation  solution  fails  to 

provide  an  adequate  description  of  the  density  distribution. 

For  n*  h  .5,  the  second  order  perturbation  solution  over-corrects 

the  first  order  solution  and  the  density  ratio  at  the  front 

itself  deviates  significantly  from  the  correct  value.  The 

asymptotic  nature  of  the  perturbation  series  is  thereby 

demonstrated  in  this  figure.  Again  the  method  of  Porzel 

yields  a  fairly  good  description  down  to  the  asymptotic  regime 

as  -»■  H-,. 
s  CJ 

The  particle  velocity  profiles  are  shown  in  Fig.  4. 

He  can  see  that  as  a  consequence  of  failing  to  conserve  mass, 
the  particle  velocity  profile  deviates  significantly  from  the 
present  numerical  solution  in  the  moderate  blast  strength 
regime  {i.e.,  o'  *8).  Porzel's  method  once  again  is  quite 
adei  i'*tt  in  the  asymptotic  regime  as 

In  general  from  the  n'  vs  plot  ei.d  the  pressure, 

density  and  particle  velocity  profiles  just  described  we  Can 
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conclude  that  Porzel's  method  Is  by  far  the  superior  of 
all  the  non-similar  analytical  techniques  and  yields  an 
adequate  description  of  the  dynamics  and  flow  structure 
of  reacting  blasts  for  the  entire  regime  of  propagation, 
I.e.,  0  <  n'  1  1. 

4*  The  Asymptotic  Regime 


In  their  analysis  of  the  asymptotic  decay  of  a 
reacting  blast  wave  to  its  C-J  state.  Levin  and  Chernyi 
obtained  an  expression  for  the  shock  decay  coefficient 

B  I 

0  =  evaluating  the  three  conservation  equations 


at  the  front  C  =  1  and  then  taking  the  limit  as  H  M 
or  e  •  l-n'  »  1  -  ->■  0,  I.e., 


CO 


Q  =  -K,€^  -  KjC 


(4.1) 


where 


K.  = 


2(7+  I)  V  , 


zJsi 


\2 


Vcj 


(4.Z) 


and 


K,  = 


7+1  <^<P 


?-i 

e-*o 


(4.3) 


The  numerical  constant  j  =  0,  1,  2  for  planar,  cylindrical 
and  spherical  waves  respectively  and  =  1/M^j®.  Levin 
and  Chernyi  then  assumed  that  the  blast  wave  solution  recovers 


the  J'T-Z  similarity  solution  in  the  limit  as  e  0,  hence 


d<p 


fi=  / 
€^0 


OO 


i  ike 


r 

K 


E) 


± 

2 


By  further  assuming  that  the  particle  velocity  gradient 


approaches  infinity  as  e  **■  0  like  one  notes  that  K  -*•  ® 

2 


like  G  ^  as  c  0  and  Eq.  4.1  yields  in  the  limit  as  e  +  0 
for  j  ^  0  the  following  relationship: 

■4 


^  =  C  e  -  ( 

where  C  is  some  arbitrary  constant.  For  the  planar  case 
where  j  -  0,  one  notes  from  Eq*  4.2  that  =  0  and  as  g  0 


the  d-T*Z  solution  yields 


? = I 
e  *  0 


yt  I 


2  I"  1  as  G  0  for  the  planar  c»5e  and  Eq*  4*1 


Hence  K 

becomes  in  the  limit  as  g  *»•  0, 


0  =  -  € 


From  the  definition  of  9,  one  notes  that 


2Rs_£2.' 

d  Rs 


-2Rc 


dV' 

dRs 


=  2R 


Cf€ 
^  (J 


since  n*  1  as  g  0.  Using  the  above  equation  and  Eq,  4,5 
for  j  =  0  and  Eq.  4.4  for  j  f  0,  c.«e  gets  the  following  resu 
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for  tHe  as’lnptotlc  decay: 


€ 


n 


for  j^o  (4*7) 


€ 


C  rn 


Rg 

Rgo 


for  ji^o  (4.8) 


where  e^*  and  C  are  arbitrary 
From  Eq.  4*7  it  is  evident  that  for 
to  the  C-0  state  is  asymptotic  with 
from  Eq.  4«8t  one  notes  that  in  the 
which  Is  finite  and  given  as 


integration  constants, 
planar  waves ^  the  decay 
R  ^  «  as  e  -►  0,  However, 
limit  as  e  +  0,  Rg  R^^ 


(4.9) 


The  present  numerical  solution  confirms  the  concTuslon  of 
Levin  and  Cherny 1  that  the  decay  to  the  C-J  state  Is 
achieved  at  a  finite  radius  for  diverging  waves  (l.e., 

J  =  Is  2).  From  the  present  numerical  solution  the  value 
of  Rgj  for  various  values  of  y  and  are  shown  in  Fig.  5. 

Here  R^j  is  calculated  from: 

a)  the  present  numerical  solution 

b)  Levin,  and  Chernyi'S  asymptotic  formula 

c)  an  asymptotic  formula  newly  derived  in  the 
present  paper 

d)  values  quoted  by  Chernyl*  Korobeinikov  et  al 


for  three  chosen  values  of 


which  characterize 


the  mixture  and  for  various  values  of  y. 
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To  check  the  asymptotic  formula  of  Levin  and 
Chernyi,  the  arbitrary  constants  and  C  in 

Eq.  4.9  are  evaluated  by  matching  Eq,  4.8,  at  some 
finite  but  small  value  of  s  (typically  1-n'  .05 

for  T}'  .95),  with  both  the  present  numerical  solution 

and  also  the  analytical  solution  using  Porzel *s  method. 

With  the  constants  thus  determined.  Levin  and  Chernyi 's 
values  of  found  from  Eg.  4.9  agree  well  with  those 
from  the  exact  numerical  solution.  Also  shown  in  Fig.  5 
are  the  values  of  ^CJ  quoted  by  Chernyi,  Korobeinikov  et 
al  '  Again,  good  agreement  with  the  present  exact 
numerical  solution  is  obtained.  The  discrepancies  that 
do  exist  are  due  to  the  choice  of  R^^  Or  (or  e^).  If 
one  chooses  closer  to  zero  (i.e.j  H,  —  M-,)  the  values 
of  Rj.j  obtained  from  the  asymptotic  formulae  will  agree 
even  better  with  the  exact  values  from  the  numerical 
solution.  Hence  we  conclude  that  the  asymptotic  formula 
given  by  Levin  and  Chernyi  predicts  correctly  the  value  of 
Rgj.  However,  from  the  present  numerical  solution,  it  can 
be  observed  from  Fig.  4  that  the  particle  velocity  profile 
as  e  -»•  0  does  not  recover  the  J^T-Z  Similarity  solution. 

The  gradient  at  the  front  as  e  0  remains  finite  Instead 
of  approaching  Infinity  as  predicted  by  the  J-T-Z  solution. 

The  value  of  the  velocity  gradient  from  the  present  exact 
solution  is  found  to  approach  that  of  the  planar  case  instead-, 
i.e. , 


d(p 

df  1  =  1 

^  6  -*0 


2 

y  + 1 


(4.1 
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From  a  physical  point  of  view  this  result  is  quite  plausible 
Since  as  the  spherical  wave  expands,  curvature  effects  become 
progressively  unimportant.  At  large  sliock  radii,  diverging 
waves  must  become  essentially  planar  when  curvature  efforts 
are  negligible. 

The  fact  that  Levin  and  Chernyi's  formula  predicts 
the  same  value  for  as  the  present  numerical  solution 
even  though  the  numerical  solution  yields  a  totally  different 
value  for  the  velocity  gradient  in  the  limit  as  -»■  R^j 
requires  clarification.  Instead  of  assuming  that  the  J-T^Z 
solution  is  recovered  as  e  -*■  0  with  the  result  that  -*■  *> 

we  take  the  particle  velocity  gradient  to  approach  e-»-0 
a  finite  constant  a  as  indicated  by  the  present  numerical 


solution.  Hence  from  Eq.  4.3  we  see  tnat  as  e 

From  Eq.  4.1  and  Eq.  4.6,  we  obtain 


0. 


<#  € 


2(  Kj€  ) 


and  carrying  out  the  integration  leads  to  the  following 
asymptotic  law  instead  of  that  given  by  Levin  and  Chernyi 
(i.e.,  Eq.  4.8) j 


(4.n) 


. . . . . . . . . . . Ill . . . . . 
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Ift  the  limit  as  e  0  and  ■*  R^j  we  get 


^ 


K2  \ 

K,  ^  ^ 


(4.12) 


instead  of  Eq.  4.9  as  given  by  Levin  and  Chernyi.  In  the 
present  asymptotic  law,  we  find  first  that  for  the  planar 
casei  J  =  0  and  Eq.  4.2  yields  =  0.  Hence  Eq.  4.12 
yields  R^j  ^  «  as  e  ->  0  as  predicted  by  Levin  and  Chernyi. 
in  arriving  at  this  result  from  Our  asymptotic  law  note  that 
it  is  not  necessary  to  assume  that  the  J-T-Z  solution  is 
recovered  in  the  limit  as  e  *  0  as  is  required  in  Levin  and 
Chernyi 's  analysis*  By  expanding  Eq.  4.12,  we  obtain 


(4.13) 


lip  to  the  order  of  From  Eq.  4.13,  we  see  therefore 

that  Rgj  is  independent  on  K^,  hence  independent  of  the 
particle  velocity  gradient.  This  independence  of  Rj.j  on 

4|-  implies  that  it  is  irrelevant  whether  the  J-T-Z 

del 

e*»-0  solution  is  recovered  or  not  in  the  limit  as  e  -»-  0. 
If  we  approach  the  exponential  term  in  Levin  and  Chernyi 's 
formula  given  by  Eq.  4.9  we  see  that  it  is  identical  to 


Eq.  4.13  to  the  order  Of  e. 


This  explains  why  identical 


results  for  h^j  are  obtained  using  Eq.  4.9  and  from  the 
present  numerical  solution  even  though  the  numerical  solution 
does  not  approach  the  J-T-Z  similarity  solution  as  e  -»■  0. 

The  present  solution  indicates  that  the  particle  velocity 


igMdient  at  the  front  approaches  that  of  the  planar  case 
;C%e. ,  Eq,  4.10),  hence  =  1 .  The  values  of  Computed 
from  Eq.  4.13  using  =  1  yield  values  much  closer  to  the 
exact  Ones  than  the  Levin  and  Cherny i  formula  even  though 
the  same  value  of  is  used  to  evaluate  the  constants  in 
the  asymptotic  formulae. 

The  range  of  validity  of  the  asymptotic  formula 

ti.e.,  Eq.  4.11  with  K  =1)  can  be  seen  from  Fig.  6  where 

2 

the  shock  decay  coefficient  3  is  plotted  against  the  shock 
Strength  n*  -  range  .6  ^  q'  £  1  the 

asymptotic  law  yields  almost  identical  results  to  the  exact 
nulserical  Solution.  In  other  words,  the  decay  from  moderate 
Strengths  to  the  final  C-J  state  is  adequately  described  by 
the  asymptotic  formulae.  Hence  as  far  as  the  decay  of  the 
reacting  blast  is  concerned  the  asymptotic  formula  coupled 
With  an  analytical  solution  such  as  the  solution  using 
Porzel's  method  can  provide  a  good  description  of  the  entire  * 
decay  process  from  n'  =  0  to  n*  =  !• 

It  is  also  interesting  to  note  from  Fig.  5  that  a 
decaying  reacting  blast  approaches  the  C-J  state  in  less 
than  an  explosion  length  for  most  gaseous  explosives  with 
the  range  5  £  iv  10  and  y  1  1.3.  From  previous  experiments 
on  direct  initiation  of  CgHg-Og  mixtures,  the  critical  energy 
for  direct  initiation  is  of  the  order  of  0.3  joules  for 
sioicniometric  mixtures  at  an  initial  pressure  of  100  torr. 

The  explosive  length  for  this  case  is  of  the  order  of  1  cm. 
Hence  the  C-J  state  is  reached  when  the  blast  radius  is  Of  the 
order  of  0.7  cm.  This  makes  experimental  verification  of  the 


-  17  - 


early  time  regime  in  which  the  reacting  blast  solutions 
are  valid  extremely  difficult,  and  existing  measurements 
of  the  shock  parameters  are  mostly  limited  to  the  asymptotic 
regime. 

5.  The  Finite  Kinetic  Rate  (F-K-R)  Solution 

For  the  numerical  solution  of  the  double  front 
F-K-R  model,  the  empirical  relationships  for  the  induction 
time  are  taken  from  White 

(og(T.-[oj^[C2Hi]» -10.81  +  (5.1) 

Strehlow 

tog(T.  ■  [02>  -  10.269  4  (5.2) 

and  Strehlow  and  Engel 

(og(T.[Oj)--  10.56  +  J2^  (5.3) 

To  illustrate  the  typical  numerical  results,  the  particular 
case  of  stoichiometric  acetyl ene*oxygen  mixture  at  an  initial 
pressure  of  100  torr  is  chosen*  Fig.  7  shows  the  variation 
of  the  shock  Hach  number  with  shock  radius  for  various  values 
of  the  initiation' energy  E^.  For  a  given  initiation  energy  E^, 
the  initial  decay  is  similar  to  that  of  the  R-B-H  solution 
which  assumes  the  reaction  rates  to  be  infinite.  This  is  to 
be  expected  since  in  the  early  time  regime  the  shock  is  very 


strong,  hence  the  induction  distance  is  sufficiently  small 
for  the  discontinuity  assumption  to  be  valid.  As  the  blast 
expands  and  weakens  the  reaction  zone  progressively  lags 
behind.  From  Fig.  7  we  note  that  the  shock  decays  past  the 
C-J  state  for  the  particular  mixture  whereupon  a  quasi -steady 
state  is  reached  when  the  shock  and  reaction  front  complex 
propagate  at  a  constant  sub  C^J  velocity*  The  existence  of 
this  quasi-steady  regime  of  propagation  has  been  confirmed 
experimentally  by  SolOukhin  Bach,  Knystautas  and  Lee  ^^2)^ 

Struck  and  Brossard  The  present  solution  also 

confirms  the  existence  of  a  critical  energy  for  direct 
initiation.  From  Fig*  7  we  note  that  for  <  .001  joules, 
the  quasi-steady  regime  cannot  be  obtained.  The  shock  decays 
continuously  to  the  acoustic  limit  and  the  trailing  reaction 
front  completely  decouples  from  the  shock  and  continues  to 
propagate  as  a  spherical  flame.  This  sub-critic.al  energy 
regime  of  propagation  has  been  studied  previously  by  the 

Authors  on  laser  spark  generated  spherical  detonations. 

The  present  numerical  solution  breaks  down,  marking  the 
termination  of  the  quasi -steady  regime  after  a  finite  period 
of  time,  this  is  due  to  the  fact  that  sharp  gradients  are 
developed  in  the  flow  field  behind  the  reaction-shock  front 
complex  during  this  quasi-steady  regime.  The  development  of 
these  sharp  gradients  in  pressure,  density  and  particle  velocity 
in  the  flow  field  are  illustrated  in  Figs.  8  to  10,  It  is  to 
be  noted  that  these  sharp  gradients  are  not  a  result  of  the 
stability  of  the  numerical  scheme  used.  Rather,  it  is  a 
physical  phenomenon  that  occurs  whenever  a  shock  wave  is 
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lUpported  solely  by  chemical  reactions  that  occur  at  a 
finite  rate.  Other  numerical  studies  of  the  propagation 
of  eKOthermiC  shock  waves  alSo  demonstrate  the 

pftset  of  oscillations  and  instability  under  the  critical 
cohditions  when  the  shock  is  mainly  supported  by  the 
energy  released  due  to  combustion.  Experimentally  it  is 
found  that  the  termination  of  the  quasi -steady  regime  Is 
marked  by  the  sudden  violent  transition  to  a  multiheaded 
spherical  detonation  wave  Therefore  as  far  as 

qualitative  features  are  concerned  the  present  numerical 
solution  confirms  existing  experimental  observations  of 
Ipherical  detonation  phenomena. 

A  quantitative  comparison  between  the  present 
solution  and  experimental  results  for  the  shock  and 
reaction  front  trajectories  are  shown  in  Figs.  11  and  12. 
in  Fig.  11,  the  stoichiometric  CgHg-Og  mixture  is  diluted 
by  70%  inert  Ng.  Initiation  is  by  laser  spark  with  spark 

energy  =  0.3  joules.  This  corresponds  to  the  sub-critical 
regime  where  the  shock  and  reaction  front  progressively 
decouple  from  each  other.  In  Fig.  12,  the  mixture  is 
undiluted  and  the  same  initiation  energy  =  0.3  joules 
corresponds  to  the  critical  energy  for  direct  initiation. 
Experimentally  it  Is  observed  that  decoupling  of  the  shock 
and  reaction  front  is  followed  by  a  quasi-steady  regime  at 
the  termination  of  which  is  the  re-establishment  to  a  highly 
asymmetrical  multi headed  spherical  detonation.  Photographic 
records  of  these  two  cases  are  given  in  a  previous  paper  (22)^ 
From  Figs.  11  and  12  we  can  see  that  the  present  solution  agrees 


quite  well  with  experiments  as  far  as  shock  and  reaction 
trajectories  are  concerned. 

For  the  same  experimental  conditions  as  Fig.  11, 
the  induction  times  for  each  particle’ crossing  the  blast 
wave  at  different  initial  radii  are  measured  from  the 
experimental  records.  A  comparison  with  the  results 
from  the  present  solution  indicates  fairly  good  agreement 
{Fig.  13).  Also  shown  In  Fig.  11  are  the  results  obtained 
from  an  analytical  solution  based  on  Porzel ‘s  method  where 
It  is  assumed  that  the  energy  released  by  chemical  reaction 
has  a  negligible  Influence  on  the  hydrodynamic  flow  structure. 
The  Importance  of  the  chemical  energy  in  perturbing  the 
non-reacting  blast  flow  structure  even  In  this  case  of  702 
Inert  nitrogen  dilution  Is  self-evident. 

From  the  present  solution  it  1$  also  possible  to 
determine  the  critical  energy  for  direct  initiation  for 
a  given  mixture  at  given  initial  conditions.  The  particular 
"case  for  stoichiometric  CgHg-Og  mixtures  is  shown  in  Fig.  14. 
The  experimental  results  are  based  on  the  laser  spark  as  an 
energy  source.  The  small  range  in  experimental  data  Is  due 
to  the  limited  power  range  of  the  ruby  laser  used.  The 
theoretical  prediction  of  the  limiting  critical  energy  based 
on  the  numerical  solution  of  the  F-K-R  model  yields  values 
which  are  considerably  lower  than  the  experimental  results 
using  the  laser  spark  from  the  Q-switched  ruby  laser.  Of  all 
the  experimental  initiation  results  currently  available  those 
with  the  laser  spark  from  the  Q-switched  System,  apart  from 


being  the  most  extensively  available,  represented  the 
shortest  deposition  times  ('^  nanoseconds)  and  hence 
Were  felt  to  be  the  ones  most  appropriate  for  comparison 
In  view  of  the  theoretical  premise  of  vanishing  pulse  length. 


We  attribute  the  discrepancy  to  two  possible  causes.  The 
first  is  that  approximate  global  kinetic  data  are  used  in 
the  theoretical  calculations.  Moreover*  such  calculations 
indicate  that  the  Solution  is  very  sensitive  to  the  specific 
values  of  the  empirical  constants  used  in  the  global  induction 
time  relations  among  which  there  exist  significant  discrepancies 
depending  on  the  source  quoting  the  relationship.  The  second 
possible  source  of  disagreement  Is  that  strictly  speaking  the 
laser  spark  from  a  Q*sw1tched  ruby  system  does  not  simulate 
ideany  the  limiting  value  of  source  energy  (i.e.,  not  infinite 
power  since  pulse  duration  20  nanoseconds).  More  reliable 
conclusions  in  this  regard  have  to  await  Initiation  experiments 
with  sub-nanOsecond  laser  sparks  which  would  simulate  the 
Infinite  source  power  and  hence  the  limiting  source  energy 
Criterion  more  closely.  We  are  currently  Initiating  such 
experiments. 

6*  Prediction  of  Transverse  Wave  Spacings 

The  cellular  dimension  in  a  multiheaded  detonation 

front  Is  an  Important  characteristic  of  the  explosive  system 

and  its  initial  conditions.  A  recent  theoretical  attempt  has 
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been  made  by  Strehlow  '  'to  predict  the  cellular  dimension  or 
the  transverse  wave  spacing.  Detailed  experimental  observations 


of  the  motion  of  a  “detonation  wavelet"  in  between 

collisions  of  transverse  waves  indicates  that  it  behaves 

like  a  decaying  reactive  blast  wave.  The  minimum  strength 

of  a  detonation  wavelet  just  prior  to  the  explosive 

re^energization  by  the  collision  of  a  pair  of  transverse 

waves  is  found  to  be  of  the  order  of  —  3.  From  Fig.  7 

we  note  that  the  minimum  shock  strength  which  corresponds 

to  an  initiation  energy  equal  to  the  critical  value  is  also 

of  the  order  of  3.  For  any  initiation  energy  less  than  the 

critical  value  the  shock  decays  to  the  acoustic  limit.  Hence 

'"it  ifc— that  the  energy  released  in  the  collision  of  a 

pair  of  transverse  waves  corresponds  to  the  critical  value 

we  can  make  a  formal  analogy  and  take  the  shock  radius  when 

Mg  is  a  minimum  as  the  longitudinal  dimension  of  the  detonation 

cell.  Using  the  present  numerical  solution,  it  is  then  possible 

to  evaluate  the  dimension  of  a  detonation  cell  based  on  the 

above  reasoning.  A  comparison  between  the  theoretically 

f31) 

predicted  values  with  the  experimental  data  of  Strehlow' 
for  argon  diluted  C2H2-O2,  H2-O2,  and  C2H^-02  mixtures  are 
Shown  in  Figs.  15  to  17  respectively.  The  predictions  of  the 
transverse  wave  dimensions  made  in  this  way  based  on  the  present 
numerical  solution  of  the  F-K-R  model  are  in  rough  agreement 
with  experimental  measurements  of  Strehlow  Once  again 

the  discrepancies  can  possibly  be  attributed  to  the  global 
approximate  kinetic  rate  assumption  used  in  the  numerical 
solution  and  the  oversimplified  representation  of  the 
re-energization  of  "detonation  wavelets"  upon  mutual  transverse 


coTTislon  by  a  point  energy  source  of  infinite  power 
density.  It  has  been  pointed  Out  by  Strehiow  that  the 
Mast  wave  model  for  transverse  wave  spacing  can  at  best 
bje  used  in  a  qualitative  sense  Since  the-energy  released  by 
the  collision  of  transverse  waves  is  not  instantaneous  as 
assumed  in  the  theoretical  model,  and  also  from  experiments 
it  is  found  that  the  center  of  the  curved  detonation  wavelet 
does  not  correspond  to  the  point  where  the  transverse  waves 
collide.  However,  in  view  of  the  fair  agreement  between 
the  theoretical  results  and  experimental  data  we  conclude 
that  the  blast  wave  model  can  be  used  to  provide  an  a-priori 
prediction  of  the  transverse  wave  spacings  at  present  until 
better  theoretical  models  are  developed.  It  is  to  be  noted 
that  the  fundamental  mechanisms  for  the  propagation  of 
multiheaded  detonation  waves  remain  obscure,  and  the 
present  blast  wave  model  fails  to  throw  light  on  why  nature 
prefers  detonative  combustion  to  proceed  in  such  a  complex 
manner  although  the  model  does  predict  the  cellular  dimensions 
themselves. 

7.  Conclusions 

Using  the  present  exact  numerical  solution  as  a 
reference  for  comparison,  it  is  found  that  Porzel 's  method 
of  the  power  law  density  profile  gives  the  most  accurate 

description  of  the  motion  of  a  reacting  blast  wave  with  the 

•  « 

widest  range  of  applicability.  In  the  asymptotic  regime  as 
the  blast  approaches  the  C-J  state,  the  present  solution 


confirms  the  results  of  Levin  and  Cherny i  in  that  for 
diverging  waves,  the  C-J  state  is  reached  at  a  finite 
radius.  Although  the  asymptotic  law  given  by  Levin  and 
Chernyi  does  give  the  correct  prediction  of  the  C-J  radius, 
the  assumption  they  invoke  that  the  J-T-Z  similarity  is 
recovered  in  the  limit  e  -►  0  is  not  supported  by  the  present 
results.  Instead  the  present  solution  indicates  that  the 
solution  approaches  the  planar  case  in  the  limit  with  all  the 
flow  gradients  behind  the  front  being  finite.  To  resolve  this 
inconsistency,  a  new  asymptotic  decay  law  is  developed  based 
on  finite  flow  gradients  behind  the  shock  as  the  C-J  state 
is  reached.  The  present  asymptotic  formula  Shows  that  to  the 

Je 

order  of  e"^,  the  C-0  radius  is  independent  Of  the  flow  gradient 
and  hence  explains  why  the  formula  of  Levin  and  Chernyi  yields 
the  correct  value  of  R^j  while  differing  in  the  assumption  as 
to  the  nature  of  the  flow  gradient  used  in  the  analysis. 

Comparing  with  the  exact  numerical  solution  it  is  found  that 
the  asymptotic  decay  law  is  valid  for  a  wide  range  of  shock 
Strengths.  Coupled  with  an  analytical  solution  for  the  Initial 
early  time  propagation  regime,  the  entire  decay  of  a  reacting 
blast  can  be  described  quite  accurately.  The  finite  kinetic 
rate  solution  is  found  not  only  to  predict  all  the  essential 
qualitative  features  of  the  propagation  of  a  spherical  combustion 
wave  but  predicts  quantitatively  the  shock  and  reaction  front 
trajectories  and  the  critical  energy  regime  for  direct  initiation. 
Using  a  blast  wave  model  for  cellular  detonations  the  theoretically 
predicted  values  of  the  transverse  wave  spacings  are  found  to  be 
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in  fair  agreement  with  the  experimental  data  of  Strelilow. 
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